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ABSTRACT 

We report the results of a 10-year timing campaign on PSR J1738+0333, a 5.85- 
ms pulsar in a low-eccentricity 8.5-hour orbit with a low- mass white dwarf com- 
panion. We obtained 17376 pulse times of arrival with a stated uncertainty smaller 
than 5 /is and weighted residual rms of 1.56 /is. The large number and precision of 
these measurements allow highly significant estimates of the proper motion /i a j = 
(+7.037 ± 0.005,-1-5.073 ± 0.012) masyr^ 1 , parallax tt x = (0.68 ± 0.05) mas and a 
measurement of the apparent orbital decay, Pb = (—17.0 ± 3.1) x 10 _15 ss _1 (all 1-a 
uncertainties). The measurements of fj, a< s and ir x allow for a precise subtraction of 
the kinematic contribution to the observed orbital decay; this results in a significant 
measurement of the intrinsic orbital decay: P 6 Int = (—25.9 ± 3.2) x 10 _15 ss _1 . This 
is consistent with the orbital decay from the emission of gravitational waves predicted 
by general relativity, P b GR = —27.71^9 x 10~ 15 ss -1 , i.e., general relativity passes 
the test represented by the orbital decay of this system. This agreement introduces 
a tight upper limit on dipolar gravitational wave emission, a prediction of most al- 
ternative theories of gravity for asymmetric binary systems such as this. We use this 
limit to derive the most stringent constraints ever on a wide class of gravity theories, 
where gravity involves a scalar field contribution. When considering general scalar- 
tensor theories of gravity, our new bounds are more stringent than the best current 
solar-system limits over most of the parameter space, and constrain the matter-scalar 
coupling constant cvq to be below the 10 -5 level. For the special case of the Jordan- 
Fierz-Brans-Dicke, we obtain the one-sigma bound Oq < 2 x 10 -5 , which is within a 
factor two of the Cassini limit. We also use our limit on dipolar gravitational wave 
emission to constrain a wide class of theories of gravity which are based on a general- 
ization of Bekenstein's Tensor- Vector- Scalar gravity (TeVeS), a relativistic formulation 
of Modified Newtonian Dynamics (MOND). 

Key words: pulsars: timing — pulsars, individual: PSR J1738+0333 — gravity: 
theories — general relativity: tests 
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Pulsar J1738+0333 is one of seven recycled pulsars discov- 
ered in a high-Galactic latitude (15° < |6| < 30°, -100° < 
I < 50° Parkes m ulti-beam pulsar survey l|jacobvl 120051 ; 
IJacobv et ai]|2007l ). It has a spin period P of 5.85 ms and it 
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is located in a low-eccentricity (e < 4 x 1CP 7 ) binary system 
with an orbital period (Pb) of 8.5 hours; the projected semi- 
major axis of the pulsar orbit (a;) is 0.3434 light seconds. 
The Parkes timing determined a phase-coherent timing so- 
lution with precise spin, orbital and astrometric parameters 
that allowed the det ection of the c ompanion of the pulsar at 
optical wavelengths (Jacobv 2005). 

The results of the optical o bservations of the com- 
panion are described in detail by lAntoniadis et all |2012l . 
henceforth Paper I); it is a white dwarf (WD), and from 
its spectrum they estimate its mass M c = 0.18lto'oo7 
(1-cr, where Mq represents one Solar mass; we also define 
m c = M c /Mq). The paper also presents measurements of 
the Doppler shifts of the WD spectral lines which are used to 
estimate the systemic radial velocity 7 = (—42 ± 16) kms" 1 
and the semi-amplitude of the WD orbital velocity pro- 
jected along the line of sight K c — (171 ± 5)kms~ 1 , con- 
sistent (b ut signifi c antly more precise) with the value ob- 
tained by IJacobvl (|2005l ). K c = (181 ± 27)kms~ 1 . This 
measurement is combined with its pulsar equivalent K v = 
2nxc/Pb — 21.10336 kms -1 to estimate the mass ratio: 
q = Mp/M c = K c /Kp = 8.1 ± 0.2. Given M c , this implies 
a pulsar mass M p = 1.46ig'o5 M©. Finally, from Kepler's 
laws they derive an orbital inclination i = 32? 6 ± 1?0. 

The position pf PSR J1738+0333 makes it detectable 
with the Arecibo Observatory's 305-m William E. Gor- 
don radio telescope, which provides about 15 times the 
sensitivity of the Parkes telescope and therefore allows 
much more precise timing. Regular Arecibo observations of 
PSR J1738+0333 started in 2003. Given the mass measure- 
ments obtained from the optical/radio data, and the knowl- 
edge that e < 4 x 10~ 7 , general relativity (GR) predicts an 
orbital decay of P 6 GR = (-27.7t\i) x ICT^ss" 1 (henceforth 
fss -1 ) due to the emission of gravitational waves (GW). Our 
early simulations suggested that using Arecibo we would be 
able to measure this effect within a few years. 

The organization of the remainder of this paper is as 
follows: in Section [5] we describe in detail the observations 
and how the pulse times of arrival and the timing solution 
were derived. In Section [3] we discuss the detection of the 
orbital decay of the system. We compare it with the pre- 
diction from GR, and derive an upper limit for any excess 
GW emission. In Section [4] we use this limit to derive an 
upper limit for the emission of dipolar GWs. In the follow- 
ing sections, we discuss the implications of this result for 
alternative theories of gravity, specifically in Section [5] for 
tensor-scalar theories of gravity and in Section [5] for a class 
of tensor- vector-s calar theori e s of gr avity similar to the the- 
ory proposed by iBekensteinl (|2004n , which is a relativistic 
formulation of the Modified Newtonian Dynamics (MOND) 
proposed bv lMilgroml (|l983l ). Finally, in Section [TJ we sum- 
marize our results and highlight the prospects opened by 
continued timing of this system. 



2 RADIO TIMING OBSERVATIONS 

2.1 Observational setup and data reduction 

PSR J1738+0333 was discovered with the 64-m Parkes tele- 
scope in an observation taken in 2001. We started timing it 
regularly on 2001 September 26 using the 21-cm multi-beam 
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Figure 1. PSR J1738+0333 average pulse profile at 1400 MHz 
taken with the WAPP spectrometers. A full 512-bin rotation cycle 
is displayed. 

receiver system and the 2 x 512 x 0.5 MHz filterbank. The 
Parkes datata king and processing are described in detail in 
IJacobvl (|2005h . 

In August 2003 we started regular observations of this 
MSP using the L-wide receiver of the 305-m William E. Gor- 
don Telescope at the Arecibo Observatory and the Wide- 
band Arecibo Pulsar Processors (WAPPs) to process the 
signal. The WAPPs are general-purpose auto-correlators 
capable of processing a b and of up to 100 MHz each 
IjDowd. Sisk fc Hagenl 120001 ). First, they digitise the incom- 
ing voltages as 3-level values. The machine then calculates 
the auto-correlation function (ACF) for a specified number 
of time lags, in our case 256. These are integrated for a 
time ts of 64 /j,s, after which they are written to disk as 16- 
bit integers for later off-line analysis. The four WAPPs were 
centered at frequencies F of 1170, 1310, 1410 and 1510 MHz. 

The ACFs from the WAPPs are, later on, read from 
disk during off-line processing. First, they are Fourier trans- 
formed using Duncan Lorimer's SIGPROC routine "filter- 
bank", producing 100 MHz-wide, 256-channel power spec- 
tra. To eliminate signal distortions caused by the edges of 
the WAPP bands we set the values of the lower and higher 
16 channels to zero, all other channel values are unchanged. 
These spectra are then divided into four sub-bands. Each 
of these is dedispersed at the DM of the pulsar. The re- 
sulting time series is then folded modulo the pulsar rota- 
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tional period, using a polynomial expansion of the pulsar 
spin phase as a function of time predicted by the best ex- 
isting ephemeris. For each of these sub-bands, the pulses 
are integrated for 4 minutes as single 512-bin pulse pro- 
files. Only on days when the pulsar is consistently faint 
(the flux density changes because of scintillation caused by 
the interstellar medium) do we integrate the pulse profiles 
for 600 seconds. All of this processing is done using the 
"sig_foldpow" routine, written by one of us (IHS). Given the 
pulsar's dispersion measure (DM, a measurement of the col- 
umn density of free electrons between the Earth and the pul- 
sar) of 33.77 cm -3 pc, the dispersive smearing per channel 
(with bandwidth SF = 100MHz/256 = 0.390625 MHz) was 
toM = 8.3^s(DM/cm- 3 pc)((S_F 1 /MHz)(F/GHz)' 3 = 68.3, 
48.7, 39.0 and 31.8 /is per channel respectively. The total 
time resolution dt = vt§~+ ^dm was therefore 93.6, 80.4, 
75.0 and 71.5 /is respectively. 

Producing a large number of such profiles has many 
advantages: a) The separation in frequency allows the use 
of polynomial coefficients derived specifically for the radio 
frequency of each sub-band in such a way that the orbital 
phase of the binary is constant for the pulse arrival time at 
each frequency, not constant for a particular time. This is 
very important for systems with short orbital periods like 
PSR J1738+0333, where the difference in arrival times at 
the lower and upper edges of the band caused by dispersion 
by the ionized interstellar medium corresponds to a non- 
negligible shift in the orbital phase, b) The production of 
many TOAs in time preserves the important orbital infor- 
mation contained in them, c) By having small bandwidths 
and integration times, we minimize profile smearing due to 
imperfections in the ephemeris and DM model. Finally, d) 
we take full advantage of the power of scintillation - in some 
occasions, the signal-to-noise ratio in a 25 MHz x 4-minute 
subsection is larger than for the whole observation. 

2.2 Derivation of times of arrival 

The best pulse profile, resulting from more than 1 hour of 
data, is displayed in Fig. [T] The narrow central component 
is the cause of the excellent timing precision for this pulsar. 
This profile does not show significant changes with radio 
frequency beyond what one should expect from the change 
in time resolution of the system dt with frequency. For this 
reason, the best profile is then cross-correlated with a ll the 4- 
minu te/25 MHz pulse profiles in the Fourier domain (|Tavlorl 
1992); the phase offset that yields the best match is used 
to derive the time of arrival (TOA) of a particular pulse 
(normally that closest to the start of each sub-integration) 
measured at the local (topocentric) time frame. In total we 
obtain 17376 TOAs from Arecibo data. The previous num- 
ber of Parkes measurements is 100. 

We carry out subseq uent TOA analysis using the 
tempo2 software package dHobbs. Edwards fc Manchester! 



2006; lEdwards. Hobbs fe Manchester] I2006T ). This uses the 
clock corrections intrinsic to the observatory, the Earth 
rotation data and the observatory coordinates to convert 
the topocentric TOAs to Terrestrial Time (TT), as main- 
tained by the Bureau International des Poids et Mesures. 
From these tempo2 derives arrival times at the Solar 
System Barycentre using the position of the Observatory 
that is calculated using the DE/LE 421 Solar System 



ephemeris (Folkner. Williams fc B oggs 2008); for details see 
lEdwards. Hobbs fc Manchester] (|2006l V These, like our tim- 
ing parameters, are expressed in Barycentric Coordinate 
Time (TCB). 

With the correct rotation count, tempo2 can vary the 
timing parameters in order to minimize the difference be- 
tween the observed and predicted barycentric TOAs (the 
timing residuals). The residuals are weighted according to 
the estimated uncertainty of each TOA. The resulting set of 
timing parameters constitutes a preliminary phase-coherent 
timing solution (an "ephemeris"). If the TOA uncertainties 
are under-estimated, or if there are systematics present in 
the data, the reduced % 2 of the fit will be larger than 1. For 
this preliminary ephemeris, the reduced x 2 is 2.05, which, 
as we show below, is mostly caused by variations in the in- 
tervening electron column density. 

The main advantage of having WAPP search data 
stored on disk and doing the analysis off-line is that it al- 
lows for iterative improvement of the pulsar ephemeris; this 
process is particularly important at the earlier stages when 
the ephemeris is not yet very precise. With each improved 
ephemeris, we dedispersed and folded the data again, obtain- 
ing pulse profiles with improved signal-to-noise ratio; which 
are then used to derive better TOAs that further improve 
the ephemeris. This avoids orbital phase-dependent smear- 
ing of the pulse profiles and the creation of orbital phase- 
dependent timing artefacts, which can corrupt the determi- 
nation of orbital parameters, particularly the orbital pha se 
and orbital period variation (N ice. Stairs fc Kasianll2008T l. 

At the early stages of our timing programme, we 
checked the timing accuracy of the WAPPs by making a 
few simultaneous observations with the Arecibo Signal Pro- 
cessor (ASP). The ASP is a real-time coherent dedispersion 
system implemente d in software, w hich can process 64 MHz 
of baseband signal (|Demoresdl2"007r i and is known to provide 
very stable timing. No significant differences in the timing 
of the two back-ends was found, this implies that, within 
measurement precision, the WAPP timing is accurate. 



2.3 DM variations and the timing model 

The preliminary ephemeris described above assumes a con- 
stant DM. As we can see in the right plot of Fig. [2] this 
ephemeris describes the TOAs rather poorly, with large 
trends in the post-fit residuals. These are especially notice- 
able once we calculate daily residual averages. 

We used this preliminary ephemeris to measure the pul- 
sar's DM for each day of observations. This is only possible 
because of the wide band of the L-wide receiver at Arecibo 
— from 1120 to 1730 MHz, of which we used the cleanest 
parts, 1120-1220 MHz and 1260-1560 MHz. The results of 
these measurements are presented in the left plot of Fig. [2] 
The DM varies significantly with time, and its variation cor- 
relates with the daily residual averages; this implies that the 
latter must, to a large extent, be caused by the former. They 
have amplitudes of the order of 0.002 cm -3 pc, which cause 
extra dispersive delays of the order of 4 /is, and a dominant 
timescale of the order of a few years. This is similar to the 
timescales associated with GWs from super massive black 
hole binaries. This highlights the importance of accurate 
and dense multi-frequency TOA measurements and precise 
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DM offset (10 cm pc) Phase ( 1 

Figure 2. Left: Daily DM measurements versus time. Given the non-linear variation of the DM, it is clear that it must be modeled 
using a high-order polynomial. We display in green the eighth-degree polynomial in Table [T] this fits the DM evolution well inside the 
range where there are DM measurements, but diverges fast outside that range. Right: Post-fit residuals versus time. The residuals of the 
Parkes (green) and Arecibo (gray) TOAs were obtained with a preliminary ephemeris that assumes a constant DM. The averages of the 
residuals for each sidereal day are indicated in red. Notice their similarity with the daily DM averages; this implies that implying the 
former are (mostly) caused by the latter. 
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Figure 3. Left: Post-fit residuals versus time. The residuals of the Parkes (green) and Arecibo (gray) TOAs were obtained with the 
timing model presented in Table[T] The averages of the residuals for each sidereal day are indicated in red. Right: Pre-fit residuals versus 
time obtained from the Westerbork Synthesis Radio Telescope TOAs at 1380 MHz (black) and 345 MHz (red). The latter residuals are 
displayed with an offset of 60 /is for clarity. These were obtained using the timing model presented in Table [l] The lack of trends in the 
Westerbork residuals implies that our timing solution describes these TOAs correctly. These data were not used in the derivation of the 
timing solution presented in this paper; their inclusion does not change any parameter significantly. 
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DM correction s when attempting; t o detect such waves with 
pulsar timing (I Verbiest et al.l 12009). 

These DM variations cannot be described by a simple 
linear variation with time. For this reason, we included eight 
DM derivatives in our definitive timing model, which is pre- 
sented in Table [T] the uncertainties were derived by tempo2 
except where stated otherwise. This DM model describes the 
daily DM measurements very well (see left plot of Fig. [2]); 
but is only valid within the time interval for which we have 
significant multi-frequency data, i.e., since the start of the 
Arecibo observations in 2003. The predicted values start di- 
verging fast outside this range, but this does not affect the 
results discussed in this paper. This method has the benefit 
that any possible correlations of the coefficients of the DM 
model with the timing parameters are taken into account, 
leading to more conservative (and in our opinion more real- 
istic) estimates of the parameter uncertainties. 

This ephemeris describes the TOAs much better than 
the preliminary ephemeris, as displayed in the left plot of 
Fig. [3] However, there are still systematic variations in the 
daily residual averages, these excursions reach a maximum 
of about 0.5 fjs and have timescales of a few months to a 
year. In order to estimate the timing parameters with more 
realistic uncertainties, we added an uncertainty of 0.5 fis in 
quadrature to the estimated uncertainty of every TOA. 

This rescaling of our errors results in a reduced \ 2 of 
1 for short timescales and 1.02 for the whole data set, sug- 
gesting that the uncertainties on the derived parameters are 
reliable. However, given the systematic nature of these ex- 
cursions, small systematics on parameters with compara- 
ble timescales (in particular the parallax) may remain. On 
the other hand, there should be no correlation with orbital 
phase, hence orbital parameters should not be affected at 
all. 

These excursions are likely to be caused by variations 
of the DM at timescales of a few months that cannot be 
taken into account by the 8-polynomial model for the DM 
variations. The alternative explanations are not as plau- 
sible: An instability in the rotation of the pulsar would 
produce excursions with larger timescales than those ob- 
served. The second derivative of the spin frequency [u = 
(—0.6 ± 2.3) x 10~ 28 Hzs -2 ] is consistent with zero, which 
suggests good long-term stability; the same is true for the 
third frequency derivative. Furthermore, the agreement be- 
tween nearby daily residual averages suggests that the tim- 
ing system used is stable. 

To further verify the integrity of this timing solu- 
tion, we compared its predictions with TOAs taken with 
the Westerbork Synthesis Radio Telescope (WSRT) in 
the Netherlands, which uses the Puma II coherent dedis- 
persi on back-end ijKaruppusamv. Stappers. fc van Stratenl 
2008). The 345 MHz data are especially useful because, in 
case of an inaccuracy in our DM model they would show 
significant trends. It is remarkable that, despite the fact 
that the model represents mostly an interpolation during the 
2009-2010 gap, no trends are discernible in either dataset 
(see right plot of Fig. [3j . 

2.4 Orbital model 

The orbit of PSR J1738+0333 has a very low eccentric - 
ity, so we use the "ELL1" orbital model (|Lange et al.ll20o"lh 



to model itQ This yields Keplerian and post-Keplerian pa- 
rameters that are very weakly correlated with each other. 
Note that, in order to estimate the "real" eccentricity of 
the binary we assumed that the Shapiro delay is as pre- 
dicted by general relativity for the values of M c and sini 
derived in Paper I. This assumption is safe because GR is 
known to provide a sufficiently accurate description of the 
distor tion of space-time around wea kly self-gravitating ob- 
jects (jBertotti. less fc Tortor^ 120031). 

According to lFreire fc Wexl (|2010P l. the orthometric am- 
plitude of the Shapiro delay (which quantifies the time am- 
plitude of the measurable part of the Shapiro delay) is, for 
this system, given by /13 = 22 ns. Fitting for this quantity 
we obtain hz = 9 ± 13 ns. This is l-cr consistent with the 
prediction but the low relative precision of this measure- 
ment implies that we cannot determine M c and sini inde- 
pendently from the existing timing data. A precise measure- 
ment of the component masses of this system from Shapiro 
delay would require an improvement in timing precision that 
is much beyond our current capabilities. 

In the right plot of Fig. 3] we display the residuals as 
a function of the orbital phase. No trends are noticeable, 
either in the residuals or their averages, this implies that the 
orbital model is not obviously flawed. This also suggests that 
the timing system is inherently stable. Furthermore, we see 
no DM variations as a function of the orbital phase (left plot 
of Fig. U); i.e., there are no obvious spurious DM artefacts 
caused by incorrect folding nor detectable dispersive delays 
in the data. 



3 RESULTS 

In what follows, we discuss some of the results of the timing 
program. In Section \3. II we briefly discuss the measurement 
of the parallax, which requires special care given the sys- 
tematics highlighted in Section \2. 31 Then in Section T3.2I we 
focus on the main result of the timing: the detection of the 
orbital decay of the system, iV We compare it with the GR 
prediction in Section 13731 



3.1 Parallax 

As mentioned in Section [231 the quoted uncertainty for the 
parallax is likely to be too small given the systematic ef- 
fects caused by uncorrected short-term DM variations. It is 
therefore important to gain a sense of whether this parallax 
estimate is accurate or if there are inconsistencies with other 
distance estimators. 

The distance we obtain from this parallax (d = 1.47 ± 
0.10 kpc) is consistent with the 1.4 kpc p redicted by the 
NE20 01 electron model of the Galaxy l|Cordes fc Laziol 
2001) for the pulsar's Galactic coordinates and DM. How- 
ever, the distance estimates based on this model have been 



1 The ELL1 timing model as implemented in the TEMPO soft- 
ware package (http:/ /sourceforge . net/projects/tempo/ ) is a mod- 
ification of the DD timing model llDamour fc Deruell jl985lll986ft 
adapted to low-eccentricity binary pulsars. In terms of post- 
Keplerian observables, it contains all those which are numerically 
relevant for systems with e <C 1 . The "Einstein delay" term is not 
relevant for such systems and is therefore not taken into account. 
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Figure 4. Left: Measured DM versus orbital phase. To the observed DMs we fit a DM variation model (green lines, see SectionfA]) using 
DM and ADM as the free parameters, and display the resulting fits for the the nominal and ±2-<r values of ADM. No significant variation 
as a function of orbital phase is detectable; indicating the lack of spurious artefacts in the data reduction and no outgassing from the 
companion. Right: Post-fit residuals versus orbital phase, which for this very low-eccentricity system is measured from ascending node 
(i.e., the mean anomaly is equal to the orbital longitude). The residuals of the Parkes (green) and Arecibo (gray) TOAs were obtained 
with the timing model presented in Table [T] In red we indicate bin averages where the bin width is 0.01 P^. No significant trends can be 
identified in the residuals or their binned averages; indicating that the orbital model can describe the orbital modulation of the TOAs 
correctly. No dispersive delays or unnacounted Shapiro delay signatures are detectable near orbital phase 0.25, nor artefacts due to 
incorrect dedispersion and folding of the data. 
© 0000 RAS, MNRAS 000, 000-000 
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Timing parameters 

Reference Time (MJD) 54600.0001776275 

Right Ascension, a (J2000) 17 h 38 m 53?9658386(7) 

Declination, <5 (J2000) 03° 33' 10'.'86667(3) 

Proper Motion (mas yr" 1 ) +7.037(5) 

Proper Motion in 8, fis (mas yr" 1 ) +5.073(12) 

Parallax, ir x (mas) (a) 0.68(5) 

Spin Frequency, v (Hz) 170.93736991146392(3) 

First Derivative of j/, v (fHz s" 1 ) -0.704774(4) 

Orbital Period P b (days) 0.3547907398724(13) 

Projected Semi-Major Axis, x (lt-s) 0.343429130(17) 

Time of Ascending Node, T asc (MJD) 54600.20040012(5) 

tjEesinu (-1.4 ± 1.1) X 10" 7 

KEecosw (3.1 ± 1.1) X 10" 7 

First Derivative of P b , P b (fss" 1 ) -17.0(3.1) 

"range" parameter of Shapiro delay, r (/is) (b) 0.8915 

"shape" parameter of Sapiro delay, s = s'mi (b) 0.53877 

Dispersion Measure, DM (cm -3 pc) 33.77312(4) 

First Derivative of DM, DM1 (cm" 3 pc yr" 1 ) +0.00108(3) 

... DM2 (cm" 3 pc yr" 2 ) -0.00041(3) 

... DM3 (cm" 3 pc yr" 3 ) -0.000279(9) 

... DM4 (cm" 3 pc yr~ 4 ) +0.000059(6) 

... DM5 (cm" 3 pc yr- 5 ) +0.0000230(19) 

... DM6 (cm" 3 pc yr" 6 ) -0.0000019(3) 

... DM7 (cm" 3 pc yr" 7 ) -0.00000090(11) 

... DM8 (cm" 3 pc yr" 8 ) -0.000000060(10) 

Test parameters 

First Derivative of x, x (fss -1 ) 0.7(5) 

Second Derivative of u, v (10 -28 Hz s~ 2 ) -0.6(2.3) 

Optical Parameters 

Companion Mass, M c (Mq) 0.181+oq^ 

Spectroscopic Companion Radius, R c (Rq) 0.037^{j 'Jjjj* 

Semi-amplitude of orbital velocity of companion, K c (kms" 1 ) 171(5) 

Derived Parameters 

Galactic Longitude, I 27? 7213 

Galactic Latitude, b 17?7422 

Distance, d (kpc) 1.47(10) 

Total Proper Motion, fj, (mas yr" 1 ) 8.675(8) 

Position angle of proper motion, M (J2000) 53?72(7) 

Position angle of proper motion, ©,j (Galactic) 116? 12(7) 

Spin Period, P (s) 0.005850095859775683(5) 

First Derivative of Spin Period, P (10" 20 ss- 1 ) 2.411991(14) 

Intrinsic P, P Int (lO-^ss" 1 ) (a) 2.243(13) 

Characteristic Age, t c (Gyr) 4.1 

Transverse magnetic field at the poles, Bq (10 9 G) 0.37 

Rate or rotational energy loss, E (10 33 erg s~ 1 ) 4.4 

Mass Function, / (M ) 0.0003455012(11) 

Mass ratio, q = M p /M c 8.1(2) 

Orbital inclination, i (°) 32.6(1.0) 

Pulsar Mass, M p (Mq) i-^to'.OS 

Total Mass of Binary, M t (M ) l-65tooe 

Eccentricity, e (3.4 ± 1.1) X 10" 7 

Apparent P b due to Shklovskii effect, P^hk ( fss -i) ( a ) 8.2l^ 

Apparent P b due to Galactic acceleration, p^al (f ss -!) ( a ) . 0.58^g'J| 

Intrinsic P b , P Int (fss" 1 ) (a) -25.9(3.2) 

Predicted P b , P 6 GR (fss" 1 ) -27.7±\% 

"Excess" orbital decay, P* s = P Int - P fc GR (fss" 1 ) (a) +2.0t 3 ^ 

Time until coalescence, r m (Gyr) ~ 13.2 



Table 1. Parameters for the PSR J1738+0333 system. In parentheses we present the 1-cr uncertainties in the last digit quoted, as 
estimated by TEMP02. If the value and uncertainty are signaled with an (a) then they were derived from a Monte-Carlo procedure 
(Section I3.4| l. (b) The Shapiro delay parameters r and s were not fitted in the derivation of the timing model; the values used were 
derived from a combination of other timing and optical parameters (Section l2.4H . All timing parameters are derived using TEMP02 and 
are displayed as measured at the Solar System Barycenter, in barycentric coordinate time (TCB). The "test parameters" were not fitted 
when deriving the main timing model, but their values were derived fitting for all the other parameters in the model. 
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shown, in some cases, to significantly under- or overestimate 
real distances, therefore this agreement with the DM pre- 
diction cannot be used as solid evidence that our parallax is 
accurate. 

This distance, when combined with the known temper- 
ature and photometric properties of the white dwarf, pro- 
duces an estimate for its radius that is very similar to that 
derived from its the spectrum (see Paper I). This suggests 
that our value for the distance is likely to be accurate given 
the present timing uncertainties. This is important because 
the distance (and the proper motion, also presented in Ta- 
ble [1]) are necessary for a correct estimate of the intrinsic 
orbital decay of the system P^"*, as discussed below. 



Following 



the 



ana lysis by 
IVerbiest. Lorimer. fc McLaughlin! (|2010l i. we find that 
there are no significant biases affecting this parallax 
measurement. 



3.2 Intrinsic orbital decay 

The intrinsic orbital decay of the system can be obtained 
from the observed orbita l period variation (Pb) by subtract- 
ing th e kinematic effects ()Shklovskiill970l : lDamour fc Tavlorl 



n Int = P b 



'->Shk 



(1) 



The same equation applies to any quantity with the dimen- 
sion of time, like the spin parameters (P, P and P Int ). 

The first term, P* cc = arPb/c is caused by the differ- 
ence of accelerations of the PSR J1738+0333 system and the 
Solar System projected along the line of sight to the pulsar, 
ax- This term is dominated by the difference of accelerations 
of the two systems caused by the average Galactic field, 
ac'- For the pulsar's Galactic coordinates of I = 27?7213 
an d b = 17? 7422 and distance, we obtain [using eq. (5) 
m iNice fc Tavlorl l| 19951 ). in combination with eq. (17) in 
Lazaridis et al. 2009] 

^Pt =0.58l^fss-\ (2) 
c 

The ax term also contains a contribution from nearby 
masses, as- Damour & Taylor (1991) present a statis- 
tical estimate of the magnitude of this effect for PSR 
B1913+16, which is dominated by large molecular clouds. 
For PSR J1738+0333, the same as would yield an orbital 
period derivative of 



as [ 



P b < 0.2 fss~ 



(3) 



which is negligible. The a s of PSR J1738+0333 is likely to 
be even smaller: although the pulsar is at a similar Galacto- 
centric distance as PSR B1913+16, it is at a Galactic height 
of ~ 0.45 kpc, which is larger than that of PSR B1913+16: 
~ 0.3 kpc. 

The second kinematic effect in eq. {TJ , commonly known 
as the "Shklovskii " effect, caused b y the centrifugal acceler- 
ation, is given by (Shklovskii 1970): 



(/'•<:, 



2 

Ma 



d 



(4) 



where fi a and fis are the proper motion in right ascension 
and declination (see Table [IJ. For the intrinsic Pb, we thus 
obtain 



A Int = -25.9 ±3.2 fss -1 . (5) 



3.3 Excess orbital decay 

From the values for q and m c in Paper I we can estimate the 
orbital decay caused by the emission of quadrupolar GWs 
for a low-eccentricity system, as predicted by GR: 

192 7T 



SGR 



(mT m c ) 5/3 



(q + 1) 1 / 3 



= -27.7tt.gfss" 



(6) 



where T fl = GM fl c" 3 = 4.925490947 its 
l|Lorimer fc Kramerl 120051 ). Subtracting this from P b Int 
(eq. ([5]l) we obtain the "excess" orbital decay relative to 
the prediction of GR, 



PF = 2.0 



+3.7 



fs S 



(7) 



This is consistent with zero. As discussed in Section [4] this 
implies that GR passes the test posed by the orbital decay of 
PSR J1738+0333. We illustrate this match in Fig. where 
we see that the mass/inclination constraints, derived from 
i& nt using eq. ((6j (i.e., assuming that GR is the correct the- 
ory of gravity), are consistent with the theory-independent0 
constraints derived from the optical observations. 



3.4 Rigorous uncertainty estimates for orbital 
decay 

To make reliable estimates of the uncertainties of these de- 
rived quantities (P 6 Gal , P b Shk , P 6 Int and P 6 XS ), we implemented 
a new Monte Carlo routine in TEMP02. From our list of 
17476 TOAs, we created 220 000 similar data sets of fake 
TOAs that have random distributions consistent with the 
original TOAs and their uncertainties. For each fake TOA 
data-set, we run tempo2 and record the resulting best-fit 
parameters to disk. We then use the parallax of each sim- 
ulation to estimate ao and calculate the derived quantities 
for that simulation using the equations above. Finally, in 
order to estimate PF, we use a random (m c ,q) pair from 
the Monte Carlo simulation in Paper I, calculate the corre- 
sponding P b GR and subtract this from P 6 Int . This procedure 
is warranted by the intrinsic lack of correlation between the 
optical and radio measurements. The computer then calcu- 
lates averages, standard deviations, medians and ±l-er per- 
centiles (presented in Table QJ for the resulting distributions 
of derived quantities. The averages and standard deviations 
are close to the estimated medians and l-cr percentiles, im- 
plying that the resulting distributions are generally close to 
Gaussian. 

With this method, we are able to take into account any 
underlying correlations between the observables, thus esti- 
mating more reliable values and uncertainties for the de- 



2 By theory-independent, in the context of this paper, we denote 
quantities which are either based on weak-fi eld gravity , which is 
known to be described extremely well by GR (Will 2006.) , or quan- 
tities which are free of any explicit strong field deviations of grav- 
ity from GR, at least within the wide class of Lorentz-invariant 
gravity theories. The mass ratio in a binary pulsar system is an 
example of the latter toamourjl2007h . 
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0.2 0.4 0.6 0.8 1 1 2 3 4 

cos i Pulsar Mass (M ) 

Figure 5. Constraints on system masses and orbital inclination from radio and optical measurements of PSR J1738+0333 and its WD 
companion. The mass ratio q and the companion mass m c are theory-independent (indicated in black), but the constraints from the 
measured intrinsic orbital decay (P^ nt , in orange) are calculated assuming that GR is the correct theory of gravity. All curves intersect, 
meaning that GR passes this important test. Left: cos i— m c plane. The gray region is excluded by the condition m v > 0. Right: m p —m c 
plane. The gray region is excluded by the condition sini < 1. Each triplet of curves corresponds to the most likely value and standard 
deviations of the respective parameters. 



rived quantities that depend only on the measured TOAs 
and their uncertainties. 



4 GENERIC TESTS OF GRAVITY THEORIES 

In order to understand the significance of the small value of 
Pb S in eq. ([7]) — the main experimental result of this paper 
— we now discuss what physical effects could in principle be 
contributing to it. According to iDamour fc Tavlorl l|l99ll ) 



P b M + Pt +PF + P b G , 



(8) 



where P^ 1 is due to mass loss from the binary, P 6 T is a con- 
tribution from tidal effects, P^ is the orbital decay caused 
mainly by the emission of dipolar GWs (and any extra 
multipole modifying the general relativistic prediction) and 
Pjp is a contribution from possible (yet undetected) varia- 
tions of Newton's gravitational constant (as measured by a 
Cavendish experiment). The first two terms are the "classi- 
cal" terms, the last two would only be non-zero for theories 
of gravity other than general relativity. 



4.1 Classical terms 

4-1.1 Mass loss 

In Appendix [X] we derive an upper limit for the mass loss 
from the companion as a function of the total mass of the 
system. For the pulsar, t he mass loss is dominat ed by the 
loss of rotational energy (|Damour fc Tavlorl Il99lh : 



M t 



E 



= 1.5 x 10~ 21 s~\ 



(9) 



which is of the same order as the upper limit for ^ . 

The contribution to the orbital varia tion due to the to- 
tal m ass loss M = M c + M v is given by |Damour fc Tavlorl 

Haas): 



pr=2^n< 0.2ft s -\ 



(10) 



which is about 20 times smaller than the current uncertainty 



4-1.2 Tidal orbital decay 

We now calculate the o rbital decay caused b y tid es. From 
eqs. (3.15) and (3.19) in lSmarr fc Blandfordl (U976D . we de- 
rive the following expression for P^: 



S7vq(q + 1) 



R c Pb sin i \ 1 



(11) 



Unlike the expressions in eq. (3.19), this equation is exact 
because it relates the synchronisation timescale r s (which 
describes the change in the companion angular velocity 
Q c , t s = — Q c /fl c ) to the timescale associated with the 
change in the orbital period (r p = Pb/Pb) assuming only 
conservation of the angular momentum. In this expression 
k = J c /(Af c Pc), where I c is the WD moment of inertia. 
White dwarfs (particularly those with a mass much below 
the Chandrasekhar limit) are sustained by the degeneracy 
pressure of non-relativistic electrons and can be well approx- 
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Figure 6. Limits on G/G and rerj derived from the measurements 
of P 6 XS of PSR J1738+0333 and PSR J0437-4715. The inner blue 
contour level includes 68.3% and the outer contour level 95.4% 
of all probability. At the origin of coordinates, general relativity 
is well within the inner contour and close to the peak of proba- 
bility density. The gray band includes regions consistent with the 
measured value and 1-cr uncertainty of G/G from Lunar Laser 
Ranging (LLR). Generally only the upper half of the diagram 
has physical meaning, as the radiation of dipolar GWs must nec- 
essarily make the system lose orbital energy. 

imated by a poly tropic sphe re with n — 1.5. For such stars, 
we have k = 0.2 ( |Motdll952h . 

The only unknown parameters in this expression are 
Q c and r s . If r 3 is much smaller than the characteristic age 
of the pulsar r c = 4.1 Gyr, then the WD rotation is already 
synchronised with the orbit (fi c = nt) and there are no tidal 
effects at all. If, on the other hand, t s > r c , then Q c can be 
much larger, but it must still be smaller than the break-up 
angular velocity D. c < sjGM c /B% = 0.038 rads" 1 . These 
conditions for Q. c and r s yield P^ < 1.4fss — 1 . Thus, even 
if the WD were rotating near break-up velocity, P^ would 
still be smaller than the uncertainty in the measurement 
of P£ s . We note, however, that the progenitor of the WD 
was very likely synchronised with the orbit. This implies 
that, when the WD formed, its rotational frequency was 
within one order of magnitude of the orbital frequency, i.e., 
s2 c < 2 x 10~ 3 rads~* (for the reasoning, see, e.g., Appendix 
B2.2 of lBassa et all ||2006f )). 



4.2 Test of GR and generic tests of alternative 
gravity theories 

The smallness of the classical terms implies that the mea- 
surement of P£ s (eq. (0) is a direct test of GR. Unlike 
many alternative theories of gravity, GR predicts P b G = 0, 
PfP = and therefore P£ B = 0. As discussed in Sec- 
tion 13.31 this is consistent with observations, which means 
that GR passes the test posed by the measurement of P£ s for 
PSR J1738+0333. In this respect PSR J1738+0333 consti- 
tutes a verification of GR's quadrupole formula with a pre- 
cision of about 15% at the 1-cr level. I n view of more strin- 
gent tests with other binary pulsars ( Kr amer et al.l 120061 ; 
IWeisberg. Nice fc Tavlo3l201Ch this result by itself does not 



seem particularly interesting. However, the large difference 
in the compactness of the two components of this binary 
system makes PSR J1738+0333 a remarkable laboratory for 
alternative gravity theories, in particular those which predict 
the emission of dipolar gravitational radiation. In Sections 
[S] and El we will confront our observations with two specific 
classes of gravity theories. In the present section, we follow a 
more generic approach, valid for gravity theories where non- 
perturbative strong-field effects are absent and higher-order 
contributions in powers of the gravitational binding energies 
of the bodies can be neglected, at least to a point where one 
does not care about multiplicative factors < 2. As an exam- 
ple, the well known Jordan-Fierz-Brans-Dicke scalar-tensor 
theory falls into this group. 

Under the assumptions above, we can write for the 
change in the orbital period caused by dipolar gravitational 
radiation damping in a low-eccentricity binary pulsar system 

P 6 D ~ -2tt n b T e m c -1_ KD S 2 + O (s 3 , c ) , (12) 

where S = s p — s c is the body-dependent term which is given 
by the difference in th e "sensitivit ies" of the pulsar, s p , and 
the companion, s c [see I Willi |l993) for the definition of s p . c ]. 
The quantity is a body-independent constant that quan- 
tifies the dipolar self-gravity contribution, and takes differ- 
ent values for different theories of gravity0 For the purpose 
of this section, we have neglected higher-order corrections 
in powers of the sensitivities in the equation above (They 
actually vanish in the Jordan-Fierz-Brans-Dicke case). The 
full non-linearity will be taken into account in Sections [5] 
and [5] (anticipating on the notation defined there, the terms 
oc Sp :C are negligible when the absolute value of the non- 
linear matter-scalar coupling constant, |/?o|, is significantly 
smaller than 2). 

The value of s p , c depends on the theory of gravity, the 
exact form of the equation of state and the mass of the 
pulsar; for a M p = 1.4 Mq neutron star, it is generally of 
the order of 0.15, the value we use in our calculations. For 
an asymmetric system like PSR J1738+0333, the sensitivity 
of the companion WD has a negligible value: in the post- 
Newtonian limit, it is given by e/M c c 2 ~ 10 ~ 4 , where t 
is the gravitational binding energy of the WD (Will 2006). 
Therefore, S = s p — s c ~ s p ^ 0, which implies that if 
kd 0, then there must be emission of dipolar GWs, and 
an associated orbital decay according to eq. (|12fl . In a double 
neutron star system we would have s p « s c and therefore 
S ~ 0, which means that we should observe fj « even if 
k d 7^ 0. It is for this reason that, despite the low relative 
precision of the radiative test in PSR J1738+0333, it rep- 
resents such a pow erful constrain t on alternative theories 
of gr avity (see, e.g. lEardlevI Il975l ; iBhat. Bailes fc Verbiestl 
l2008h . Apart from this, the use of optical data is very im- 
portant because they provide estimates of q and m c that are 
free of explicit strong field effects — unlike in the case of the 
binar y pulsar PSR B1913+16 (see lWeisberg. Nice fc Tavlorl 
l2010h. or for many of the parameters of the double pulsar 
l|Kramer et~ai]|2006h . 

3 In general scalar-tensor theories of gravity, we have k,£> = 
V (1 - 7 PPN )~\ where r) = 4/3 PPN - 7 PPN - 3 is the Nordtvedt 
parameter, a combination of PPN parameters related to the vio- 
lation of the strong equivalence principle. 
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The last contribution to P™ comes from a possible 
contribution to the orbital change by a varying gravita- 
tional constant (P G in eq. ©). In the worst case, P° 
and P 6 G could both be large (in violation of GR) but just 
happen to cancel each other in the PSR J1738+0333 sys- 
tem because of different signs. To disentangle these effects 
there are two methods. First, one can use the best cur- 
rent limits from tests in the Solar System, notably Lunar 
Laser Rang ing (LLR), which yields G/G = (- 0.7 ± 3.8) x 
10" 13 yr _1 (|Hofmann. Miiller fc BiskupektoloD , and obtain 
for PSR J1738+0333 a (conservative) upper limit of 

Pi = -1%Pb = (+0.14 + 0.74) fss" 1 , (13) 
ijDamour. Gibbons fc Tavlorlll988l :lD amour fc Tavlorlll99ir i. 



Therefore, P, u = P h x 



-P b G = 1.9^; 8 7 fss- 



a typical sensitivity s p = 0.15 
(-0.8 ± 1.6) x 10" 4 , 



KD 



which yields, for 



(14) 



a limit that is a factor o f eight more stringent than the limit 

from PSR J1012+5307 (|Lazaridis et al.ll2009h. 

The second method, developed in ILazaridis et al.l 

(2009), combines two binary pulsar systems with different 
orbital periods. The method is based on the fact that a wide 
orbit is more sensitive to a change in the gravitational con- 
stant but less affected by the emission of dipolar GWs, in 
comparison to a more compact orbit. If we combine the P^ s 
of PSR J1738+0333 with that of a binary pulsar with a 
longer orbital period we obtain a simultaneous test for k_d 
and G. 

When calculating P G for a combined limit on kd and 
G based on two binary pulsars, we need to account for mass 
variations in compact stars as a result of a changing gravita- 
tional constant. Otherwise our limit on G will be too tight 
(Nordtvedt 1990). As a first approximation, that only ac- 
counts for the influence of the local value of G, we can use 
eq. (18) in Nordtvedt (1990): 



pG 



-2° 
G 



2q + 3 



3(7 + 2 



(15) 



2q + 2 ' 2q + 2 

As in cq. (|12[) . the contribution from the sensitivity of 
the white-dwarf companion, s c , can be neglected. For 
PSR J1738+0333, the correction factor due to the sensitiv- 
ities (i.e., the parenthesis on the right hand side of eq (|15[) ) 

is about 0.85; 

As in ILazaridis et all (I2009T), we use the P? 13 of 
PSR J0437-4715 (|Deller et al.1 |200S| ; IVerbiest et all I2008T I 
to complement our measurement (see eq. (J7J)). 

PSR J0437-4715 has a slightly higher mass than 
PSR J1738+0333, and we will account for this in the sensi- 
tivity by having s p scal e proportional to the mass, as sug - 
gested by eq. (B.3) of iDamour fc Esposito-Faresd l|l992h . 
The joint probability density function for G/G and kd is 
displayed in Fig. [5] At the origin of coordinates, GR is well 
within the inner 68% contour and close to the peak of prob- 
ability density, i.e., it is consistent with the experimental 
results from these two binaries. Marginalizing this probabil- 
ity distribution function, we obtain 



G/G 



= (-0.6 + 1.6) x 10 
= (-0.009 ± 0.022) H 
= (-0.3 + 2.0) x 10~ 4 , 



-12 -1 



(16) 
(17) 



LLR 



l«ol O 

A 




LLR 

Jl 141-6545 
J1738+0333 



Figure 7. Solar-system and binary pulsar \-a constraints on the 
matter-scalar coupling constants cto and /3o- Note that a log- 
arithmic scale is used for the vertical axis |ao|, i.e., that GR 
(«o = A) = 0) ' s sen * a * an infinite distance down this axis. 
LLR stands for lunar laser ranging, Cassini for the measure- 
ment of a Shapiro time-delay variation in the Solar System, and 
SEP for tests of the strong equivalence principle using a set 
of neutron star- white dwarf low-eccentricity binaries (see text). 
The allowed region is shaded, and it includes general relativity. 
PSR J1738+0333 is the most constraining binary pulsar, although 
the Cassini bound is still better for a finite range of quadratic cou- 
pling No- 



where Ho is Hubble's constant l|Riess et al.l [2009) and 
the uncertainties are l-cr. The P™ measurement of 
PSR J0437— 4715 is mostly responsible f or the limit on 
G/G, and it has therefore not improved since lLazaridis et alJ 
(2009). The P£ s measurement of PSR J1738+0333 is mostly 
responsible for the limit on kd, which has improved by a fac- 
tor of ~ 6 since Lazaridis et al. (2009). Although the limit 
on G/G derived from binary pulsar experiments is one order 
of magnitude less restrictive than that derived from LLR, it 
is of interest because it represents an independent test. 

The analysis presented in this section is restricted to 
gravity theories that do not develop nonperturbative strong- 
field effects in neutron stars. This assumption is well justi- 
fied for PSR J1738+0333, since such effects do not seem to 
exist in other binary pulsars with similar masses, or even 
with a higher mass like in the case of PSR J1012+5307 
l|Lazaridis et alj|2009r ). Even when non-perturbative effects 
do develop, we will show below that the higher-order correc- 
tions entering eq. (|12|l do not change the conclusions quali- 
tatively. 
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5 CONSTRAINTS ON SCALAR-TENSOR 
THEORIES OF GRAVITY 

The most natural alternatives to GR involve a scalar field ip 
contributing to the gravitational interaction, in addition to 
the metric describing usual spin-2 gravitons. In these 
theories, matter is assumed to be universally coupled to 
a physical metric g^ v = A 2 ((p)g^, v , where A(ip) is a non- 
vanishing function defining the matter-scalar coupling. 

It is convenient to expand this function around the 
background value ipo imposed by the cosmological evolution, 
as In A(ip) = lnA((po)+aa{ip-<PQ) + ?2l3o(<p-ipo) 2 + - ■ • , where 
ao defines the linear matter-scalar coupling constant, /3o the 
quadratic coupling of matter to two scalar particles, and 
we will not consider higher-order vertices in the following 
(those corresponding to «o and /3o are diagrammatically rep- 
resented on the axes of Figs. [7] and [8] the circles meaning a 
matter source ) . GR corresponds to ap = /3p = 0, and Brans - 
Dicke theory (jjordanll 1953 ; iFierjl 19561 ; iBrans fc Dickelll96ll ) 
to a 2 , = 1/(2w B d + 3) and /3 = 0. 

The predictions of such theories have been 

carefully studied in the lit eratu re (|Willl 1 19931 ; 

IDamour fc Esposito-Faresel Il992l . 1 1991 Il996al lbl Il998l ). 
In strong-field conditions, notably within and near a neu- 
tron star, the coupling constants ao and f3o are modified by 
self-gravity effects, and become body-dependent quantities, 
a a and Pa (A being a label for the body), which can be 
computed by numerical integration of the field equations. 
One needs to assume a specific equation of state (EOS) 
for nuclear matter in such integrations, and we will use 
the moderate one of |5 amour fc Esposito-Faresel i|l996bl ) 
in the following. The dependence on the stiffnes s of th e 
EOS is illustrated in IDamour fc Esposito-Faresel (|l998t ) . 
but this does not change the relative strength of the various 
binary-pulsar tests. 

The body-dependent parameters «a and /3 a enter all 
observable predictions, and for instance, the effective gravi- 
tational constant between two bodies A and B reads 

G A B = G,A 2 (ip {] )-(l + OLAa B ), (18) 

where G* denotes Newton's bare constan10. This induces a 
generic time dependence of the gravitational constant (cf. 
eq. (I15[l ). as well as a violation of the strong equivalence 
principle (SEP): The acceleration of a body depends on its 
gravitational binding energy. Let us also quote the expres- 
sions taken by the genera lizations of Eddi ngton's PPN pa- 
rameters 7 PPN and /3 PPN (|Eddingtonlll923l ). as the first will 
differ in Section [6] below: 

, CtAOLB 



JAB 



Pic 



= 1 + 



1 + CtACtB 
1 PAOtBCtC 



2(1 + «Aas)(l + ctAac) ' 



(19) 
(20) 



where A, B, C denote a priori three bodies, but B — C is 
allowed. The relativistic periastron advance, proportional to 



(2 + 2 7 * 



ft ) in the PPN formalism, becomes now a 



combination of the above expressions, explicitly written in 
eq. (9.20a) of IDamour fc Esposito-Faresel l| 19921 ). 

But the most spectacular deviation from GR is that 



4 Newton's constant G, as measured in the Cavendish experi- 
ment, is given by G,A 2 (ipo) ■ (1 + a p)- 



scalar waves are now also emitted by any binary system, 
thus contributing to the observed variation of the orbital 
period. For asymmetric systems, notably the neutron star- 
white dwarf binary studied in the present paper, the main 
contribution comes from dipolar waves: 



-2im, 



G, M c 



q l + e 2 /2 
q + 1 (1-e 2 ) 5 / 2 



n2 



(21) 



where the eccentricity e is negligible in our case. [The lowest- 
order expansion of (a p — oi c )' 2 in powers of the sensitivities 
s PiC is denoted as k,dS 2 in eq. (I12|l above. In the present sec- 
tion, we are numerically taking into account the full nonlin- 
ear dependence on the bodies' self-gravity.] The companion's 
scalar charge a c w ao because of its small binding energy, 
while the pulsar's scalar charg e a v may be of order 1 in 
some theorie s even if qq ~ (|Damour fc Esposito-Faresel 
1 19931 . fl996bl ). The orbital decay from dipolar gravitational 
wave emission (eq. (|2ip ). which is of order 0(l/c 3 ), is thus 
generically much larger than the usual quadrupole of order 
C(l/c 5 ). An observed Pt consistent with general relativity 
therefore strongly constrains scalar-tensor theories. 

This is illustrated in Fig. [7] where we also display the 
1-a constraints imposed by Solar-System tests (the Cassini 
measurement of the Shapiro delay, see Bertotti, less & Tor- 
tora 2003 and the aforementioned LLR experiment, Hof- 
mann, Miiller & Biskupek 2010) and binary pulsars, accord- 
ing to the latest literature (Weisberg, Nice & Taylor 2010 
for PSR B1913+16, Stairs et al. 2002 for PSR B1534+12, 
Kramer et al. 2006 for PSR J0737-3039A/B, Bhat, Bailes & 
Verbiest 2008 for PSR J1141-6545 and Gonzalez et al. 2011 
for the SEP test with an ensemble of binary pulsars with 
wide orbits and low eccentricities). For PSRs B1913+16 and 
B1534+12, we multiplied the error bars on their measured 
P 6 Int by two because of the known uncertainties on their dis- 
tances, which preclude an accurate estimate of the kinematic 
contributions to their orbital decay. 

Note that the LLR constraints present a (deformed) 
vertical asymptote at /?o = — 1, whereas SEP and dipolar- 
radiation dominated constraints (from PSRs J1738+0333 
and J 114 1-6545) exhibit one for /3 ~ -1.5. This differ- 
ence comes from the higher-order corrections in powers of 
the sensitivities, that we neglected in eq. (|12p above but 
which are taken into account in the present section. These 
higher-order terms are vanishingly small in the Earth-Moon 
system relevant to LLR, but they are numerically significant 
for neutron stars, and can even dominate the lowest-order 
contributions. This is notably the case when /?o = — 1 — a%, 
which implies kd = and a dipolar radiation actually start- 
ing at order 0(s 4 ) instead of 0(s 2 ). However, the proximity 
of these two vertical asymptotes illustrates that Fig. 7 would 
keep a similar shape even if we neglected all higher-order cor- 
rections. This would just slightly shift the various curves, as 
would also do a different choice of the nuclear EOS. This 
justifies a posteriori the lowest-order truncation used in our 
analysis of Section 14.21 above, even for non-negligible val- 
ues of \Po\- For large values of this coupling constant, the 
higher-order terms are responsible for the dissymmetry of 
Fig. 7 with respect to the sign of 1 + /3o- 

A comment on the companion mass used in eq. (|21[) . 
The mass of the white dwarf companion is derived from the 
optical data of Paper I, which yields GM C , and not G*M C . 
The difference is a neglibible factor I + Qqi which we anyway 
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take into account when calculating the constraints on the 
matter-scalar coupling constants ao and /?o- 

Figure [7] shows that scalar-tensor theories with a 
quadratic matter-scalar coupling /3q < —5 are forbid- 
den, whatever the value of the linear coupling ao. This is 
due to the nonperturbativ e strong-field effects studied in 
iDamour fc Esposito-Faresel (|l993l . Il996bh . For (3 Q > -5, 
the limits on ao are now derived either from the Cassini 
experiment or from PSR J1738+0333. For positive /3o 
Solar System tests used to provide the best constraints 
on qq, but this has recently changed: PSR J1141— 6545 
( Bhat . Bailes fc Verbiestl |2008) is more constraining than 
the Solar System tests for j3 > 7 and PSR J 1738+0333 is now 
the most constraining of all for /3o > 0.7. The same is true for 
the —4.8 < Po < — 2.4 range. The special case /3q = (the 
Jordan-Fierz-Brans-Dicke theory of gravity) is in the region 
where the Cassini experiment is still more sensitive. Our 1-a 
pulsar limit a 2 < 2 x 10 -5 converts into ojbd > 25000. This 
is within a factor of 1.7 of the precision of the Cassini exper- 
iment. We obtain the same constraint i n the massive Brans - 
Dicke theories recently considered in lAlsing et all (|201ll ) 
when the scalar's mass m s c 2 < h/Pt — 1.35 x 10~ 19 eV 
(where h is Planck's constant), and no longer any signifi- 
cant constraint for larger scalar masses, consistently with 
Fig. 1 of that reference. 

Overall, PSR J1738+0333 provides significantly better 
constraints than the previous best binary pulsar experiment, 
PSR J1141-6545 (Fig. If the limits obtained with that 
or other systems improve in the near future that would rep- 
resent an important confirmation of the results obtained in 
this paper. 



6 CONSTRAINTS ON TEVES-LIKE THEORIES 

A tensor-vec tor-scalar (TeVeS ) theory of gravity has been 
proposed by iBekensteinl (|2004n to account for galaxy rota- 
tion curves and weak lensing without the need for dark mat- 
ter. This is a relativistic realizatio n of the modif ied Newto- 
nian dynamics (MOND) proposal (Milgrom 1983), which in- 
troduces a fundamental acceleration scale ao ~ 10 -10 ms~ 2 
(not to be confused with the matter-scalar coupling con- 
stant ao defined above). One of the difficulties is to be able 
to predict significant deviations from Newtonian gravity at 
large distances, while being consistent with Sola r-System 
and binary-pulsar tests of GR at sma ller scales l|Sandersl 
ll997l ; lBruneton fc Esposito-Faresdl2007r ). Indeed, the scalar- 
field kinetic term of TeVeS is an unknown nonlinear func- 
tion, which must take different forms at small and large 
distances. This function can have a natural shape only if 
ao| > \J rotz/rMOND w 0.05, where tqu is the orbital radius 
of Uranus and tmond = \/GM e /a ~ 7000 AU « 0.1 lt-yr, 
otherwise the model would predict ano malies too large to b e 
consistent with planetary ephemerides l|Laskar et alJ feoog). 
Below this value, the function needs to be tuned, and even 
fine-tuned for much smaller |ao|, and it merely cannot exist 
any longer if |qo| < tqu / ?"mond ~ 0.003 (it would need to 
be bi- valued). Using binary pulsars to constrain the matter- 
scalar coupling constant ao within TeVeS is thus particularly 
interesting. 

In the following, we shall not take into account all 
the subtle details of TeVeS, which are not relevant for 
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Figure 8. Similar theory plane as in Fig. [7] but now for the (non- 
conformal) matter-scalar coupling described in the text, general- 
izing the TeVeS model. Above the upper horizontal dashed line, 
the nonlinear kinetic term of the scalar field may be a natural 
function; between the two dashed lines, this function needs to 
be tuned; and below the lower dashed line, it cannot exist any 
longer. The allowed region is shaded. It excludes general relativ- 
ity (ao = /3o = 0) because such models are built to predict modi- 
fied Newtonian dynamics (MOND) at large distances. Note that 
binary pulsars are more constraining than Solar-System tests for 
this class of models (and that the Cassini bound of Fig.[7]does not 
exist any longer here). For a generic nonzero /3o, PSR J1738+0333 
is again the most constraining binary pulsar, while for 0q rj 0, 
the magnitude of |ao| is bounded by the J0737— 3039 system. 

our conclusions. In particular, Bekenstein introduced spe- 
cific scalar- vector couplings to avoid superluminal propaga- 
tion, but this i s actu al ly not necessary to respect c ausal- 
ity llBrunetonl 120071: iBruneton fc Esposito-Faresel 120071 : 
iBabichev. Mukhanov fc Vikmanl 120081 '). We will thus focus 
on the small-distance behavior of this theory, and assume 
that the scalar-field kinetic term takes its standard form 
in this limit. We will also neglect the contributions of the 
vector field to the dynamics and the gravitational radiation, 
because they depend again on some coupling constant which 
can be chosen small enough. Let us just underline that if the 
vector field carries away some significant amount of energy, 
then binary-pulsar data are even more constraining than 
what we obtain below. Neglecting these contributions is thus 
a conservative choice. To simplify, we will here assume that 
the vector field of TeVeS is aligned with the proper time 
direction of matter. 

The crucial difference of TeVeS, with respect to the 
standard scalar-tensor theories of Section [S] is the expres- 
sion of the physical metric g M „ to which matter is universally 
coupled. In the rest frame of matter, it reads goo = A 2 ((p)goo 
for the time component, but g t j — A~ 2 (i fi)gg for the s patia l 
ones. The actual model constructed in IBekensteinl (|2004h 
assumes A(ip) = exp(ao(p), but we shall here generalize it 
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and write In A(tp) = In A((p ) + a (ip - (p ) + \Po{f — ^o) 2 , 
as in Section [5] above. We are thus considering here a 
class of TeVeS-lik e theories rather than the single model of 
IBekensteinl i|2004h . 

This specific form of the physical metric results in sev- 
eral important differences with respect to Section [S] First 
of all, to compute the matter-scalar coupling parameters 
ctA and /3a corresponding to a strongly self-gravitating 
body A, one must numerically integrate new field equa- 
tions. To s ave space, we just quote here the dif ferences with 
respect to iDamour fc Esposito-Faresel (|l996bl ). where the 
standard scalar-tensor case was discussed in detail: all fac- 
tors A 4 ((p) in its Eqs. (3.6a,b,d,e,h) and (3.14d,h) should be 
replaced by A~ 2 (ip), while the factor A 3 (ip) of its eq. (3.6f) 
should become A~ 3 (ip). Finally, the (e — 3p) term of its 
eq. (3.6d) and (3.14d) should be replaced by the sum (i+3p). 
[On the other hand, note that the (i — p) term of its 
eq. (3.6d) does not change, nor the sign of a(tp) in its 
eq. (3.6 e).] These modifications a r e con sistent with the re- 
sults of lLaskv. Sotani fc Gianniosl l|200ct ). After numerically 
integrating these field equations, starting from central values 
of the scalar field and matter density, one can extract the 
mass of the neutro n star and its effective scalar charge a a in 
the same way as in lDamour fc Esposito-Faresel (Il996bh . The 
effective quadratic coupling /3a is then defined as 9«a/cVo 
for a fixed baryonic mass. For some binary-pulsar tests, it is 
also necessary to compute the star's moment of inertia, and 
its derivative with respect to the background scalar field <po 
at constant baryonic mass. 

The form of the physical metric also changes the post- 
Newtonian predictions of the theory. Instead of eqs. (|I8|) and 
(O above, we now get 

Gab = G,(l + a A a B ), (22) 

PPN , / 00 \ 

Iab = 7 =1, (23) 

while eq. (|20|) is unchanged. Note that the effective gravita- 
tional constant Gab does not de pend any longer on A 2 (ipo). 
In the actual TeVeS theory of IBekensteinl (|2004h . where 
Po — 0, the weak-field effective gravitational constant is 
thus a true constant given by G*(l + an), in particular 
time-i ndependent as was first derived in iBekenstein fc Sagil 
(2008). This constancy also explains why LLR and other 
SEP tests do not constrain at all the theories with /3o = 0, 
as illustrated in Fig. [8] 

The expression for the periastron advance takes the 
same form as in standard scalar-tensor theories in terms 
of 7ab and Pbc, but now with ^ab = 1. In the weak- field 
conditions of the Solar System, the fact that 7 PPN = 1, like 
in GR, explains why the Cassini time-delay experiment does 
not provide any constraint on TeVeS-like theories. 

All the other predictions of scalar-tensor theories, in- 
cluding the dipolar damping given by eq. (|2I[) . keep the same 
forms in terms of the body-dependent quantities a a and Pa, 
but their numerical values do differ because of the numeri- 
cal integrations described above. For /3o = 0, i.e., the actual 
TeVeS theory of Bekenstein (2004), one finds that a a = 
cto independ ently of the body. This can be de duced from 
eq. (4.14) of IDamour fc Esposito-Faresel fl996a|), where the 
scalar field is no longer sourced by the trace of the energy- 
momentum tensor, as in standard scalar-tensor theories, but 
by the same combination (e + 3p) as the gravitational po- 
tential. The constancy of a a = ao, when /3n = 0, is also 



confirmed by our numerical integrations. The consequence 
of this extra symmetry of TeVeS is that the dipolar radia- 
tion predicted by eq. (|2ip merely vanishes, and asymmetrical 
neutron star-white dwarf binaries like PSR Jf 738+0333 no 
longer provide any stringent constraints. This is illustrated 
in Fig. [5] where several binary-pulsar tests now present a 
bump on the vertical axis, /3n = 0. As compared to Fig. [7] the 
main effect of the non-conformal matter-scalar coupling as- 
sumed in TeVeS-like models is thus to displace these bumps 
from Po ~ — 1 or —2 to the axis fio — 0. 

For generic TeVeS-like models with |/3o| > 0.1, we find 
again that PSR JI738+0333 is the most constraining binary 
pulsar. We also note that several binary pulsar tests are 
now more constraining than all Solar-System experiments 
(i.e., the LLR line of F ig. [S] while the perihelion shift of 
Mercury <|Shapirolll99(ih gives even weaker constraints, and 
all deflection or time-delay predictions coincide with those 
of GR at the first post-Newtonian order). 

From a theoretical point of view, the most natural value 
of the linear matter-scalar coupling constant |an| would be 
of order unity. Figure [S] shows this is experimentally for- 
bidden, even on the vertical axis fto = thanks to sev- 
eral binary pulsars whose tests do not rely on the magni- 
tude of the dipolar radiatiorQ. They actually impose that 
|qo| < 0.035 < ^/fQjj/rMOND, therefore TeVeS needs to as- 
sume an unnatural shape of the function defining its scalar's 
dynamics. 

A possible way t o avoid fine tuning in TeVeS-like models 
has be en proposed in lBabichev. Deffavet fc Esposito-Faresel 
where it was shown that |ao| = 1 can be consis- 
tent with Solar-System and binary-pulsar tests, thanks to a 
screening of all scalar-field effects at small distances. There- 
fore, the results of the present paper do not rule out TeVeS, 
but show that its original 2004 formulation by Bekenstein 
may need to be amended. At present, even its original writ- 
ing is consistent, although it does need some tuning. 



7 CONCLUSIONS AND PROSPECTS 

It is quite fortunate that the timing precision of 
PSR JI738+0333 allows a precise measurement of the key 
observables necessary for an estimation of the intrinsic or- 
bital decay (Pt, /J. a , Us and n x ) and that the optical obser- 
vations provide a precise estimate of a general relativistic 
prediction for the orbital decay, thus allowing a stringent 
limit on the emission rate for dipolar GWs. As a result of 
this, PSR JI738+0333 is already the most constraining bi- 
nary pulsar for (conformally-coupled) scalar-tensor theories 
of gravity. 

PSR J 1 738+0333 is also the most constraining test of 
TeVeS-like theories when the quadratic matter-scalar cou- 
pling constant /3o|>0.f. In fact, for /?o < —1 and /3o > 3, 
such theories are excluded altogether. Bekenstein's TeVeS 
(a special case with f3o = 0) is still allowed by the results 
of this experiment, but already needs some tuning given the 
small limit |ao| < 0.035 that we obtain from the double 

5 The two "horns" of the constraints imposed by PSR B1913+16 
come from the fact that this system provides only three post- 
Keplerian observables, i.e., only one test, and some fine-tuned 
models can thus pass it. 
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pulsar results ijKramer et alj [20061 ). We note that the pre- 
cision of the latter result has greatly improved since 2006 
and will be presented in a forthcoming publication (Kramer 
et al., in prep.). This will significantly reduce the allowed 
values of |ao| in the gap around /3o = 0. As a consequence, 
all surviving TeVeS-like theories will have to be unnaturally 
fine-tuned, including Bekenstein's TeVeS. 

Continued timing of PSR J1738+0333 will give us a 
much more precise measurement of Pb and (to a smaller ex- 
tent) of ir x . This will allow a further improvement of the 
estimate of the intrinsic orbital decay. However, unless we 
can make a more precise determination of m c , we will not be 
able to improve upon the current estimate of -P 6 GR ; this im- 
plies that the precision of this test is unlikely to improve by 
more than a factor of 2. Despite that, a precise measurement 
of P 6 Int is still very useful: together with the measurement 
of q, it will eventually allow (with the assumption of the 
correctness of GR) a very precise estimate of the system 
masses and a precise calibration of the mass-radius relation 
for white dwarfs. 
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Figure Al. Upper limits for mass loss from the white dwarf 
derived from upper limit of pulsar spin-down energy incident on 
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APPENDIX A: CONSTRAINING THE MASS 
LOSS OF THE WD COMPANION 

The mass loss of the companion, M c , can be constrained 
by assuming (safely) that the energy driving the process 
comes from the pulsar. Its rotational energy loss is given by 
E = -IpQpflp = 4TY 2 IpPi nt P' 3 = 4.4 x 10 33 ergs-\ where 
Pint is its intrinsic spin-down (which is calculated from the 
observed spin-down after taking into account kinematic ef- 
fects, see Section 1321) Ip ~ 10 45 gem 2 its moment of inertia 
and Q p its angular velocity. Some of this radiated energy 
impacts the white dwarf; its maximum occurs when the pul- 
sar is beaming all the power right on the orbital plane, in 
which case the white dwarf intercepts a fraction given by 
/ = 2R c /(2na), where R c is the WD companion radius, 
0.037 R© (see Paper I) and a — xc(q+ l)/sini is the orbital 
separation. This then provides the escaping particles with a 
velocity vo such that, through conservation of energy 

\m c v 2 q <E^^M c < 2.60 x 10~ 23 s" 1 ( ~ 2 M t , (Al) 
2 na V c / 
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where we express the mass loss rate as a fraction of the 
total mass of the system. To first approximation, this initial 
velocity has to be larger than the escape velocity from the 
surface of the WD, v e = ^2GM C /R C = 1.37 x 10 3 kms" 1 = 
0.0046 c; this immediately yields M c /M t < 1.2 x 10~ 18 s _1 . 

In what follows, we assume that the gas escapes from 
the WD equally in all directions. The matter density is then 
given by p(r) = M c /[^nr 2 v(r)\ and the electron density by 



n(r) 



Mc 



4:TTr 2 X[ipv(r) ' 



(A2) 



where r is the distance to the center of the WD, fip is the 
proton mass and X is the number of nucleons per electron 
(1 for Hydrogen, 2 for heavy elements). Therefore, for any 
particular density profile n(r), the larger the plasma velocity, 
the larger is M c . As calculated below, for the latter to be 
of any practical significance in the present case, vo has to 
be at least of the order of a few percent of the speed of 
light c. In this case the plasma would be little affected by 
the WD escape velocity and have a nearly constant velocity. 
This results in an upper limit for the density profile given 
by n(r) = M c / (Aitr 2 X [ipva). Therefore, the upper limit for 
the extra column density along the length I of the line of 
sight from the pulsar to Earth is given by 



(A3) 



H = / n(r)dl=— : 

J-r p coai, 4nX^ P v rpsmtp 

where r v is the distance between the pulsar and the WD, 
and ip is the angle between the direction to the pulsar and 
the direction to the observer as seen from the WD, in such 
a way that r = \{r v sin-i/Q 2 + l 2 ] 1 ^ 2 . This e quatio n is similar 
to eq. (5) of iKaspi. Tauris fc Manchester! (|l996l V 

For a very low-eccentricity orbit r p ~ a. For superior 
and inferior conjunction ip — n/2 + i and ip = n/2 — i, 
respectively, therefore 



A_ff = H sup — Hint = tt~ 



M c 



2ty Xfipvo asini 



(A4) 



Dividing by the total mass of the binary Mt = M p + M c we 
obtain 



M t 



Up xccot i . „ 
M c i 



= 1.97 x 10" 



- (?) 



ADM. 



(A5) 



In this expression 



use AH = ADM x 3.0857 x 



10 ls cm pc 1 to convert the column density to a DM, X = 2 
as an upper limit for X and the x and i measured for this sys- 
tem (see Tabled]). The DM variations for PSR J1738+0333 
are displayed in the left plot of Fig. [4] Fitting an equivalent 
of eq. (|A3|) to these DMs (also displayed in the same plot) 



we obtain ADM = (0.24 ± 1.09) x 10~ 4 cm" 3 pc, which is 
basically consistent with no observed dispersive delays. 

Using the upper limit represented by eq. (|A1|) and the 
1-ct upper limit (84% C. L.) for ADM in eq. ![A"5"|) we obtain 
^ < 2.5 x 10 -21 s _1 , which happens at a velocity of about 
0.1 c (see Fig.[M|. This is large enough for the velocity to be 
nearly constant as a function of r (eq. IA5[) but small enough 
for the Newtonian expression for the kinetic energy (eq. IA1|) 
to provide a good approximation. 

It is possible that the gas is being ejected slowly 
compared to the pulsar wind, in which case it might be 



blown away from the pulsar by its wind, forming a highly 
anysotropic "cometary" tail. In that case the low orbital 
inclination of this system would mean that the plasma em- 
anating from the WD might avoid crossing the line of sight, 
producing no extra dispersive delays detectable at the Earth. 
In such a case we can only rely on the energy balance to con- 
strain the mass loss, which at these low plasma velocities 
allows, in the extreme case discussed above, for very sub- 
stantial mass loss. There are, however, other binary pulsars 
with companion WDs of similar mass li ke PSR J0751+1807 
jNice et al.ll2005l ) and PSR J1909-3744 jjacobv et al.ll2005h 
that happen to be observed at much higher inclinations, and 
in no case is such a cometary tail ever observed. 
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